# appendix.R
# Yimeng Li
# This script produces Figures A3-A6 in
#   "Relaxing the No Liars Assumption in List Experiment Analyses".
# First version: Dec 20, 2018. This version: Dec 20, 2018.

# Set working directory:
# setwd("")

# Load libraries:
library(bindata) # required by sim_bounds.R
library(dplyr)
library(ggplot2)
library(nleqslv) # required by list_relaxed_liars.R
library(tidyr)

source("sim_bounds.R")

#*********************************************************
# More Simulations: High vs. Low Correlation
#*********************************************************

# Simulation 2.1: Add 0.05 correlation to every pair of items to Simulation 1.1

probvec.2.1 <- c(1/6, 1/2, 2/3, 2/3, 0.62)

corrmat.2.1 <- cbind(c(1,0.05,0.05,0.05,0.05),
                     c(0.05,1,0.05,0.05,0.05),
                     c(0.05,0.05,1,0.05,0.05),
                     c(0.05,0.05,0.05,1,0.05),
                     c(0.05,0.05,0.05,0.05,1))

bounds.2.1.no <- sim_bounds(probvec.2.1, corrmat.2.1, J=4, affirmative = TRUE, seed=2000, liardist="no liars", listtype="Low Corr.")
bounds.2.1.const <- sim_bounds(probvec.2.1, corrmat.2.1, J=4, affirmative = TRUE, seed=2001, liardist="constant liars", listtype="Low Corr.")
bounds.2.1.nonconst <- sim_bounds(probvec.2.1, corrmat.2.1, J=4, affirmative = TRUE, seed=2002, liardist="nonconstant liars", listtype="Low Corr.")

# Simulation 2.2: Add 0.05 correlation to every pair of items to Simulation 1.2

probvec.2.2 <- c(0.5, 0.5, 0.15, 0.85, 0.62)

corrmat.2.2 <- cbind(c(1,-0.55,0.05,0.05,0.05),
                     c(-0.55,1,0.05,0.05,0.05),
                     c(0.05,0.05,1,0.05,0.05),
                     c(0.05,0.05,0.05,1,0.05),
                     c(0.05,0.05,0.05,0.05,1))

bounds.2.2.no <- sim_bounds(probvec.2.2, corrmat.2.2, J=4, affirmative = TRUE, seed=2003, liardist="no liars", listtype="Low Corr. Des.")
bounds.2.2.const <- sim_bounds(probvec.2.2, corrmat.2.2, J=4, affirmative = TRUE, seed=2004, liardist="constant liars", listtype="Low Corr. Des.")
bounds.2.2.nonconst <- sim_bounds(probvec.2.2, corrmat.2.2, J=4, affirmative = TRUE, seed=2005, liardist="nonconstant liars", listtype="Low Corr. Des.")

# Simulation 2.3: Add 0.1 correlation to every pair of items to Simulation 1.1

probvec.2.3 <- c(1/6, 1/2, 2/3, 2/3, 0.62)

corrmat.2.3 <- cbind(c(1,0.1,0.1,0.1,0.1),
                     c(0.1,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.2.3.no <- sim_bounds(probvec.2.3, corrmat.2.3, J=4, affirmative = TRUE, seed=2006, liardist="no liars", listtype="Mid Corr.")
bounds.2.3.const <- sim_bounds(probvec.2.3, corrmat.2.3, J=4, affirmative = TRUE, seed=2007, liardist="constant liars", listtype="Mid Corr.")
bounds.2.3.nonconst <- sim_bounds(probvec.2.3, corrmat.2.3, J=4, affirmative = TRUE, seed=2008, liardist="nonconstant liars", listtype="Mid Corr.")

# Simulation 2.4: Add 0.1 correlation to every pair of items to Simulation 1.2

probvec.2.4 <- c(0.5, 0.5, 0.15, 0.85, 0.62)

corrmat.2.4 <- cbind(c(1,-0.5,0.1,0.1,0.1),
                     c(-0.5,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.2.4.no <- sim_bounds(probvec.2.4, corrmat.2.4, J=4, affirmative = TRUE, seed=2009, liardist="no liars", listtype="Mid Corr. Des.")
bounds.2.4.const <- sim_bounds(probvec.2.4, corrmat.2.4, J=4, affirmative = TRUE, seed=2010, liardist="constant liars", listtype="Mid Corr. Des.")
bounds.2.4.nonconst <- sim_bounds(probvec.2.4, corrmat.2.4, J=4, affirmative = TRUE, seed=2011, liardist="nonconstant liars", listtype="Mid Corr. Des.")

# Simulation 2.5: Add 0.15 correlation to every pair of items to Simulation 1.1

probvec.2.5 <- c(1/6, 1/2, 2/3, 2/3, 0.62)

corrmat.2.5 <- cbind(c(1,0.15,0.15,0.15,0.15),
                     c(0.15,1,0.15,0.15,0.15),
                     c(0.15,0.15,1,0.15,0.15),
                     c(0.15,0.15,0.15,1,0.15),
                     c(0.15,0.15,0.15,0.15,1))

bounds.2.5.no <- sim_bounds(probvec.2.5, corrmat.2.5, J=4, affirmative = TRUE, seed=2006, liardist="no liars", listtype="High Corr.")
bounds.2.5.const <- sim_bounds(probvec.2.5, corrmat.2.5, J=4, affirmative = TRUE, seed=2007, liardist="constant liars", listtype="High Corr.")
bounds.2.5.nonconst <- sim_bounds(probvec.2.5, corrmat.2.5, J=4, affirmative = TRUE, seed=2008, liardist="nonconstant liars", listtype="High Corr.")

# Simulation 2.6: Add 0.15 correlation to every pair of items to Simulation 1.2

probvec.2.6 <- c(0.5, 0.5, 0.15, 0.85, 0.62)

corrmat.2.6 <- cbind(c(1,-0.45,0.15,0.15,0.15),
                     c(-0.45,1,0.15,0.15,0.15),
                     c(0.15,0.15,1,0.15,0.15),
                     c(0.15,0.15,0.15,1,0.15),
                     c(0.15,0.15,0.15,0.15,1))

bounds.2.6.no <- sim_bounds(probvec.2.6, corrmat.2.6, J=4, affirmative = TRUE, seed=2009, liardist="no liars", listtype="High Corr. Des.")
bounds.2.6.const <- sim_bounds(probvec.2.6, corrmat.2.6, J=4, affirmative = TRUE, seed=2010, liardist="constant liars", listtype="High Corr. Des.")
bounds.2.6.nonconst <- sim_bounds(probvec.2.6, corrmat.2.6, J=4, affirmative = TRUE, seed=2011, liardist="nonconstant liars", listtype="High Corr. Des.")

bounds_A3 <-
  rbind(bounds.2.1.no, bounds.2.1.const, bounds.2.1.nonconst,
        bounds.2.3.no, bounds.2.3.const, bounds.2.3.nonconst,
        bounds.2.5.no, bounds.2.5.const, bounds.2.5.nonconst) %>%
  mutate(liardist = factor(liardist, levels = c("nonconstant liars", "constant liars", "no liars")),
         listtype = factor(listtype, levels = c("Low Corr.", "Mid Corr.", "High Corr."))) %>%
  gather(key="bound_type", value="bound_value", -c("listtype", "liardist"))

FigureA3 <- ggplot(data=bounds_A3, aes(x=liardist, y=bound_value, fill = bound_type)) +
  geom_violin(scale = "area", position=position_dodge(0.5)) +
  stat_summary(fun.y=median, na.rm = TRUE, geom="point", size=2, color="black") +
  geom_hline(yintercept = 0.62, linetype = "dashed") +
  facet_grid(listtype ~ ., scales="free") +
  labs(x = "Liar Distribution", y = "Bounds Estimates", fill = "") +
  scale_y_continuous(limits = c(-0.1, 1.1), breaks = c(0, 0.25, 0.5, 0.62, 0.75, 1)) +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(text = element_text(size=12)) +
  coord_flip() +
  scale_fill_manual(values=c("#CCCCCC", "#999999"))

ggsave("./FigureA3.png", FigureA3, width = 8, height = 4.8, units = "in", dpi = 600)

bounds_A4 <-
  rbind(bounds.2.2.no, bounds.2.2.const, bounds.2.2.nonconst,
        bounds.2.4.no, bounds.2.4.const, bounds.2.4.nonconst,
        bounds.2.6.no, bounds.2.6.const, bounds.2.6.nonconst) %>%
  mutate(liardist = factor(liardist, levels = c("nonconstant liars", "constant liars", "no liars")),
         listtype = factor(listtype, levels = c("Low Corr. Des.", "Mid Corr. Des.", "High Corr. Des."))) %>%
  gather(key="bound_type", value="bound_value", -c("listtype", "liardist"))

FigureA4 <- ggplot(data=bounds_A4, aes(x=liardist, y=bound_value, fill = bound_type)) +
  geom_violin(scale = "area", position=position_dodge(0.5)) +
  stat_summary(fun.y=median, na.rm = TRUE, geom="point", size=2, color="black") +
  geom_hline(yintercept = 0.62, linetype = "dashed") +
  facet_grid(listtype ~ ., scales="free") +
  labs(x = "Liar Distribution", y = "Bounds Estimates", fill = "") +
  scale_y_continuous(limits = c(-0.1, 1.1), breaks = c(0, 0.25, 0.5, 0.62, 0.75, 1)) +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(text = element_text(size=12)) +
  coord_flip() +
  scale_fill_manual(values=c("#CCCCCC", "#999999"))

ggsave("./FigureA4.png", FigureA4, width = 8, height = 4.8, units = "in", dpi = 600)

#*********************************************************
# More Simulations: High vs. Low Prevalence
#*********************************************************

#------------------------Prevalence = 0.62------------------------------

# Simulation 3.1: Same as Simulation 1.3

probvec.3.1 <- c(1/6, 1/2, 2/3, 2/3, 0.62)

corrmat.3.1 <- cbind(c(1,0.1,0.1,0.1,0.1),
                     c(0.1,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.3.1.no <- sim_bounds(probvec.3.1, corrmat.3.1, J=4, affirmative = TRUE, seed=3000, liardist="no liars", listtype="High Prev.")
bounds.3.1.const <- sim_bounds(probvec.3.1, corrmat.3.1, J=4, affirmative = TRUE, seed=3001, liardist="constant liars", listtype="High Prev.")
bounds.3.1.nonconst <- sim_bounds(probvec.3.1, corrmat.3.1, J=4, affirmative = TRUE, seed=3002, liardist="nonconstant liars", listtype="High Prev.")

# Simulation 3.2: Same as Simulation 1.4

probvec.3.2 <- c(0.5, 0.5, 0.15, 0.85, 0.62)

corrmat.3.2 <- cbind(c(1,-0.5,0.1,0.1,0.1),
                     c(-0.5,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.3.2.no <- sim_bounds(probvec.3.2, corrmat.3.2, J=4, affirmative = TRUE, seed=3003, liardist="no liars", listtype="High Prev. Des.")
bounds.3.2.const <- sim_bounds(probvec.3.2, corrmat.3.2, J=4, affirmative = TRUE, seed=3004, liardist="constant liars", listtype="High Prev. Des.")
bounds.3.2.nonconst <- sim_bounds(probvec.3.2, corrmat.3.2, J=4, affirmative = TRUE, seed=3005, liardist="nonconstant liars", listtype="High Prev. Des.")

#------------------------Prevalence = 0.27------------------------------

# Simulation 3.3: Prevalence of the sensitive item = 0.27, otherwise identical to Simulation 1.3

probvec.3.3 <- c(1/6, 1/2, 2/3, 2/3, 0.27)

corrmat.3.3 <- cbind(c(1,0.1,0.1,0.1,0.1),
                     c(0.1,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.3.3.no <- sim_bounds(probvec.3.3, corrmat.3.3, J=4, affirmative = TRUE, seed=3006, liardist="no liars", listtype="Low Prev.")
bounds.3.3.const <- sim_bounds(probvec.3.3, corrmat.3.3, J=4, affirmative = TRUE, seed=3007, liardist="constant liars", listtype="Low Prev.")
bounds.3.3.nonconst <- sim_bounds(probvec.3.3, corrmat.3.3, J=4, affirmative = TRUE, seed=3008, liardist="nonconstant liars", listtype="Low Prev.")

# Simulation 3.4: Prevalence of the sensitive item = 0.27, otherwise identical to Simulation 1.4

probvec.3.4 <- c(0.5, 0.5, 0.15, 0.85, 0.27)

corrmat.3.4 <- cbind(c(1,-0.5,0.1,0.1,0.1),
                     c(-0.5,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.3.4.no <- sim_bounds(probvec.3.4, corrmat.3.4, J=4, affirmative = TRUE, seed=3009, liardist="no liars", listtype="Low Prev. Des.")
bounds.3.4.const <- sim_bounds(probvec.3.4, corrmat.3.4, J=4, affirmative = TRUE, seed=3010, liardist="constant liars", listtype="Low Prev. Des.")
bounds.3.4.nonconst <- sim_bounds(probvec.3.4, corrmat.3.4, J=4, affirmative = TRUE, seed=3011, liardist="nonconstant liars", listtype="Low Prev. Des.")

bounds_A5 <-
  rbind(bounds.3.1.no, bounds.3.1.const, bounds.3.1.nonconst,
        bounds.3.2.no, bounds.3.2.const, bounds.3.2.nonconst,
        bounds.3.3.no, bounds.3.3.const, bounds.3.3.nonconst,
        bounds.3.4.no, bounds.3.4.const, bounds.3.4.nonconst) %>%
  mutate(liardist = factor(liardist, levels = c("nonconstant liars", "constant liars", "no liars")),
         listtype = factor(listtype, levels = c("High Prev.", "High Prev. Des.", "Low Prev.", "Low Prev. Des."))) %>%
  gather(key="bound_type", value="bound_value", -c("listtype", "liardist"))

FigureA5 <- ggplot(data=bounds_A5, aes(x=liardist, y=bound_value, fill = bound_type)) +
  geom_violin(scale = "area", position=position_dodge(0.5)) +
  stat_summary(fun.y=median, na.rm = TRUE, geom="point", size=2, color="black") +
  geom_hline(yintercept = 0.27, linetype = "dashed") +
  geom_hline(yintercept = 0.62, linetype = "dashed") +
  facet_grid(listtype ~ ., scales="free") +
  labs(x = "Liar Distribution", y = "Bounds Estimates", fill = "") +
  scale_y_continuous(limits = c(-0.1, 1.1), breaks = c(0, 0.27, 0.5, 0.62, 0.75, 1)) +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(text = element_text(size=12)) +
  coord_flip() +
  scale_fill_manual(values=c("#CCCCCC", "#999999"))

ggsave("./FigureA5.png", FigureA5, width = 8, height = 6.4, units = "in", dpi = 600)

#********************************************************************
# More Simulations: Affirmative vs. Negative Sensitive Responses
#********************************************************************

#------------------------Affirmative Sensitive Response------------------------------

# Simulation 4.1: Same as Simulation 1.3

probvec.4.1 <- c(1/6, 1/2, 2/3, 2/3, 0.62)

corrmat.4.1 <- cbind(c(1,0.1,0.1,0.1,0.1),
                     c(0.1,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.4.1.no <- sim_bounds(probvec.4.1, corrmat.4.1, J=4, affirmative = TRUE, seed=4000, liardist="no liars", listtype="Aff. Sens.")
bounds.4.1.const <- sim_bounds(probvec.4.1, corrmat.4.1, J=4, affirmative = TRUE, seed=4001, liardist="constant liars", listtype="Aff. Sens.")
bounds.4.1.nonconst <- sim_bounds(probvec.4.1, corrmat.4.1, J=4, affirmative = TRUE, seed=4002, liardist="nonconstant liars", listtype="Aff. Sens.")

# Simulation 4.2: Same as Simulation 1.4

probvec.4.2 <- c(0.5, 0.5, 0.15, 0.85, 0.62)

corrmat.4.2 <- cbind(c(1,-0.5,0.1,0.1,0.1),
                     c(-0.5,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.4.2.no <- sim_bounds(probvec.4.2, corrmat.4.2, J=4, affirmative = TRUE, seed=4003, liardist="no liars", listtype="Aff. Sens. Des.")
bounds.4.2.const <- sim_bounds(probvec.4.2, corrmat.4.2, J=4, affirmative = TRUE, seed=4004, liardist="constant liars", listtype="Aff. Sens. Des.")
bounds.4.2.nonconst <- sim_bounds(probvec.4.2, corrmat.4.2, J=4, affirmative = TRUE, seed=4005, liardist="nonconstant liars", listtype="Aff. Sens. Des.")

#------------------------Negative Sensitive Response------------------------------

# Simulation 4.3: Sensitive response is negative, otherwise identical to Simulation 1.3

probvec.4.3 <- c(1/6, 1/2, 2/3, 2/3, 0.38)

corrmat.4.3 <- cbind(c(1,0.1,0.1,0.1,0.1),
                     c(0.1,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.4.3.no <- sim_bounds(probvec.4.3, corrmat.4.3, J=4, affirmative = FALSE, seed=4006, liardist="no liars", listtype="Neg. Sens.")
bounds.4.3.const <- sim_bounds(probvec.4.3, corrmat.4.3, J=4, affirmative = FALSE, seed=4007, liardist="constant liars", listtype="Neg. Sens.")
bounds.4.3.nonconst <- sim_bounds(probvec.4.3, corrmat.4.3, J=4, affirmative = FALSE, seed=4008, liardist="nonconstant liars", listtype="Neg. Sens.")

# Simulation 4.4: Sensitive response is negative, otherwise identical to Simulation 1.4

probvec.4.4 <- c(0.5, 0.5, 0.15, 0.85, 0.38)

corrmat.4.4 <- cbind(c(1,-0.5,0.1,0.1,0.1),
                     c(-0.5,1,0.1,0.1,0.1),
                     c(0.1,0.1,1,0.1,0.1),
                     c(0.1,0.1,0.1,1,0.1),
                     c(0.1,0.1,0.1,0.1,1))

bounds.4.4.no <- sim_bounds(probvec.4.4, corrmat.4.4, J=4, affirmative = FALSE, seed=4009, liardist="no liars", listtype="Neg. Sens. Des.")
bounds.4.4.const <- sim_bounds(probvec.4.4, corrmat.4.4, J=4, affirmative = FALSE, seed=4010, liardist="constant liars", listtype="Neg. Sens. Des.")
bounds.4.4.nonconst <- sim_bounds(probvec.4.4, corrmat.4.4, J=4, affirmative = FALSE, seed=4011, liardist="nonconstant liars", listtype="Neg. Sens. Des.")

bounds_A6 <-
  rbind(bounds.4.1.no, bounds.4.1.const, bounds.4.1.nonconst,
        bounds.4.2.no, bounds.4.2.const, bounds.4.2.nonconst,
        bounds.4.3.no, bounds.4.3.const, bounds.4.3.nonconst,
        bounds.4.4.no, bounds.4.4.const, bounds.4.4.nonconst) %>%
  mutate(liardist = factor(liardist, levels = c("nonconstant liars", "constant liars", "no liars")),
         listtype = factor(listtype, levels = c("Aff. Sens.", "Aff. Sens. Des.", "Neg. Sens.", "Neg. Sens. Des."))) %>%
  gather(key="bound_type", value="bound_value", -c("listtype", "liardist"))

FigureA6 <- ggplot(data=bounds_A6, aes(x=liardist, y=bound_value, fill = bound_type)) +
  geom_violin(scale = "area", position=position_dodge(0.5)) +
  stat_summary(fun.y=median, na.rm = TRUE, geom="point", size=2, color="black") +
  geom_hline(yintercept = 0.62, linetype = "dashed") +
  facet_grid(listtype ~ ., scales="free") +
  labs(x = "Liar Distribution", y = "Bounds Estimates", fill = "") +
  scale_y_continuous(limits = c(-0.1, 1.1), breaks = c(0, 0.25, 0.5, 0.62, 0.75, 1)) +
  theme_bw() +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(text = element_text(size=12)) +
  coord_flip() +
  scale_fill_manual(values=c("#CCCCCC", "#999999"))

ggsave("./FigureA6.png", FigureA6, width = 8, height = 6.4, units = "in", dpi = 600)
